A numerical method for determining the interface free energy 
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We propose a general method (based on the Wang-Landau algorithm) to compute numerically 
free energies that are obtained from the logarithm of the ratio of suitable partition functions. As 
an application, we determine with high accuracy the order-order interface tension of the four-state 
Potts model in three dimensions on cubic lattices of linear extension up to L = 56. The infinite 
volume interface tension is then extracted at each j3 from a fit of the finite volume interface tension 
to a known universal behavior. A comparison of the order-order and order-disorder interface tension 
at /3c provides a clear numerical evidence of perfect wetting. 
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I. INTRODUCTION 

Interfaces play a major role in various physical phe- 
nomena in, e.g.. Statistical Mechanics, Soft Condensed 
Matter, Particle Physics and Biology. Of particular inter- 
est are fine (or rough) interfaces, whose long range fluc- 
tuations are determined by massless modes. For those 
systems, the infrared properties are universal and can 
be described by general models like the capillary wave 
model [r or the Nambu-Goto string [21 [3]. 

An interesting class of interfaces are those related to 
the free energy of topological excitations. This is the 
case, e.g., for the order-order interface in the Potts model 
or for the tension of the 't Hooft loop in SU(A^) gauge 
theories. The free energy of the topological object for 
a system on a (hyper) cubic domain of linear size L can 
be extracted from the ratio of the partition function of 
the system in the presence of the topological excitation 
(which can be enforced using suitable boundary condi- 
tions) over the partition function of the system with pe- 
riodic boundary conditions. In more detail, if Z{L) is 
the partition function in the presence of a topological ex- 
citation and Z{L) the partition function with periodic 
boundary conditions, Fi , the free energy of the interface 
(which we assume to be translationally invariant in one 
direction) is given by 



F/(i) = -log 



Z{L) 
Z{L) 



log(L) 



The interface tension a is then obtained as 
,. Fj{L) 



(1) 



(2) 



with D the dimensionality of the system. 

With some noticeable exceptions, it is not known how 
to determine from first principles the analytical behav- 
ior of interfaces as a function of the couplings of the sys- 
tem. Moreover, it is a notoriously hard problem to access 
directly partition functions in Monte Carlo simulations, 
since these quantities have exponential fluctuations in the 
volume (see e.g. [HIS])- Most of the solutions adopted in 



the literature (e.g. [SHH]) consist in relating the inter- 
face tension (or its derivative) to quantities that can be 
reliably determined via Monte Carlo simulations. How- 
ever, these methods generally introduce large systematic 
and/or statistical errors. Hence, a direct determination 
of the interface free energy as ratio of partition functions 
is desirable. 



II. THE MODEL 

The partition function Z of a system at a tempera- 
ture T = 1/(^^/3), with ks the Boltzmann constant, is 
obtained as the integral (or the sum for discrete levels) 
over the energy E of the density of states g{E) weighted 
with the Boltzmann factor e~^^: 



Z^ f g{E)e-^^dE 



(3) 



Monte Carlo methods sample efficiently the distribution 
g{E)e~^^ , and are best suited for determining statisti- 
cal averages of observables with Gaussian fluctuations. 
An independent strategy for studying statistical proper- 
ties of a system consists in the numerical determination 
of g{E). This can be achieved using the Wang-Landau 
algorithm [3]. In this article, we propose a method to 
extract the interface tension using Eqs. (IT]) and (pi) based 
on the Wang-Landau algorithm. The method is tested 
on the four-state Potts model in three dimensions. 

In the q-state Potts model the fundamental degrees of 
freedom arc spin variables that can take the integer values 
0, . . . , q — 1. The Hamiltonian of the model computed on 
a configuration q is given by 






(4) 



where J is the strength of the interaction, S^.g. is the 
Kronecker delta function of the spin variables qi,qj on 
neighbor sites i,j and the sum {ij) is over nearest neigh- 
bors. For a system of finite size, periodic boundary condi- 
tions in all directions are imposed. The partition function 



is then given by 






-PE 



(5) 



where the first sum is over all possible configurations q 
and the second over all allowed energies (from now on, 
we redefine /3 as J/{kBT)). At zero temperature, there 
are q stable vacua . In two and three spatial dimensions, 
the system transitions from the low temperature ordered 
phase, in which the spins are predominantly in one of the 
q values, to the high temperature disordered phase at the 
critical temperature Tc- 

For simplicity, we now specialize to the three- 
dimensional case. At zero temperature, it is possible to 
enforce an interface separating two regions with two dif- 
ferent vacua by imposing twisted boundary conditions 
in one direction (e.g. the third direction). We consider 
here only twists of one unit, i.e. twists for which spins 
qi at points i with coordinates {xi,X2,L) are replaced 
by {qj + 1) mod q, where j = (a;i,X2,0) . We call H{q) 
the corresponding Hamiltonian. The configuration that 
minimizes the energy has a misalignement of the spins 
by one unit on a plane orthogonal to the third direction. 
At finite temperature, this rigid interface between the 
two vacua can fiuctuate and near the phase transition 
becomes dominated by massless modes (rough phase). 
The partition function for the system with an interface 
is given by 



y^g-/9ff(9) 



J29iE)e~ 



-liE 



(6) 



where the tilde indicates that those quantities have to 
be computed for the system with Hamiltonian H. With 
these definitions, for a system on a cubic lattice of size L, 
the free energy and the tension of the interface between 
two ordered states (order-order interface) are given re- 
spectively by Eq. ^ and Eq. ^. 



III. THE NUMERICAL DENSITY OF STATES 

To access directly the partition functions ^ and (61), 
we use the Wang-Landau algorithm ^,. This algorithm 
modifies directly the density of states by performing a 
random walk in energy space. A random update of a 
spin is accepted with a probability mhY{l,g{E)/g[E')}, 
where E and E' are respectively the energies before and 
after the update. After the update, g{E) is modified, s.t. 
g{E) — )■ kg{E), where now E is the energy of the configu- 
ration after the update. The update satisfies the detailed 
balance in the limit A; — > 1. However, starting from a 
k ^ 1 is important to obtain a first rough approximation 
of g{E). For this reason, we begin the simulation with 
a relatively large value of k, namely fc = 3. Then, when 
g{E) has converged, we reduce k — >■ yk and repeat the 
cycle until k is small enough for the systematic errors to 



be significantly smaller than the statistical errors. Each 
cycle defines one iteration of the algorithm. We study 
cubic lattices of size L = 16, 20, 24, 28, 32, 40, 48, 56. For 
each volume, we perform 20 independent simulations for 
both periodic and twisted boundary conditions. For 
the smallest lattice, we start from a constant density of 
states. For the other lattices, we use as an input an in- 
terpolation of the density of states determined on the 
system with the closest smaller size. After the simula- 
tion, g{E) is normalised so that f/(_E,„in) = log(4) and 

9{En,in) = log(4L). 

In a Wang-Landau-type simulation, the convergence of 
the density of the states and the saturation of the error 
present potential issues |10H14j . The original implemen- 
tation [9] used a flatness criterion for the histogram of the 
visits to the various energy levels: when the histogram is 
fiat within some tolerance, we assume that g{E) at that 
level of iteration has converged. However, the tolerance 
is somewhat arbitrary. In ^10] . it was proposed that the 
density of state converges when each energy value is vis- 
ited at least l/y/\ogk times. Increasing the number of 
visits will not decrease the statistical error. We will refer 
to this proposal as the Zhou-Bhatt convergence criterion. 
This criterion however does not address the convergence 
of the measured density of state to the true density of 
state. In [131 HI! i* was shown that possible systematic 
errors due to the convergence to the wrong density of 
states can be eliminated if we require a number of visits 
at least equal to 1/logfc for each energy level (Morozov- 
Lin convergence criterion). 

In Fig. [l] we provide a comparison of the fiat histogram 
criterion, the criterion proposed in llOj and the criterion 
introduced in [T3] for the specific heat 



C = 



t 
L3 



{E^)-{Ey) 



(7) 



at /3 = 0.316 on a 16'^ lattice. The figure shows that af- 
ter 20 iterations all the three criteria are at convergence. 
Moreover, the three estimates for the specific heat are the 
same within errors, the typical size of the errors being a 
few percent. This general feature is independent from /3. 
An explicit comparison of the Zhou-Bhatt criterion with 
the Morozov-Lin criterion focusing again on the specific 
heat is provided in Fig. [2]for a wide range of /?. This sug- 
gests that for observables that are accurate to the level 
of the percent, among which are interface tensions, the 
three criteria yield compatible results. In terms of statis- 
tical errors, the Zhou-Bhatt criterion seems to yield the 
largest error bars. However, the Zhou-Bhatt algorithm 
converges to the density of the states in a CPU time that 
is a factor of 16 smaller than the original Wang-Landau 
and a factor of 64 smaller than the Morozov-Lin crite- 
rion. This study suggests that for our application the 
Zhou-Bhat criterion (i.e. number of visits ex l/\/log k) is 
adequate from the numerical point of view for the level 
of precision requested by our study. Hence, we used this 
criterion to decide when a given iteration had converged. 
Based on the study of Fig.[l] we performed 23 iterations. 




Figure 1: (Color online) Comparison between three different criteria to assess the convergence of the density of states, for the 
specific heat Cv at /3 = 0.316. The horizontal line is the central value obtained with the Morozov-Lin criterion (middle plot) at 
the 27th iteration. The left plot has been obtained with the Wang-Landau flat histogram method, while the right plot shows 
the data obtained with the Zhou-Bhatt criterion. The vertical line marks the 23rd iteration. 
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Figure 2: (Color online) The difference between the specific 
heat determined with the Zhou-Bhatt {Cy) and the Morozov- 
Lin criterion {Cy) divided by Cy (both taken after 23 itera- 
tions), on a wide range of /3 for the 16'^ lattice. 



which is a conservative estimate of the number of itera- 
tions needed for the convergence of the algorithm. 

To perform further tests of our implementation, we cal- 
culated some thermodynamical quantities like the critical 
temperature, the latent heat and the entropy density, and 
compared our results with Refs. [TS1[TS]. Following p^ . 
we use three different estimators for the transition tem- 
perature on lattices with finite extension. The first, /S^t, 
is defined to be the value for which the canonical distri- 
bution P{E, /3) = g{E)e~^^ has two equal maxima. The 
second (/3^) is the position of the central energy of the 



latent heath, i.e. the value of /3 satisfying the equation 
1 



eif,)^-[e+iP') + e-{f)] , 



(8) 



where e^ (e — E/N, with N — L^) are the locations of 
the maxima of the probability distribution at 0^. The 
third (/3j?) is the location of the maximum of the specific 
heat [7] All the critical temperatures are extrapolated to 
the infinite volume limit using the ansatz 



pi{L) - /?: + 



L3 



(9) 



We perform the extrapolation by using the six largest 
lattices. The extrapolated values agree in the infinite- 
volume limit (see Fig. Isl) . Averaging over the three de- 
terminations gives Pc = 0.3143103(9). The latent heat 
per site Ae can be calculated from the maxima of the 
specific heat [17] : 



C„,ax(i) = c+J(/3,)2(Ae)2L3 



(10) 



Our estimate is Ae — 1.16454(16). Finally, the numeri- 
cally determined entropy density 



/3(e-/) 



(11) 



with / = — log Z/{f3N) the free energy density, is plotted 
in Fig. [4] All those quantities are always less than two 
standard deviations from the corresponding determina- 
tions of [15], which have been obtained on larger lattices. 



IV. INTERFACE TENSIONS 

We now move to the discussion of our results for the 
interface tension. In the ordered phase, an interface can 
form between two regions of space that are in two differ- 
ent vacua. This interface is called the order-oder inter- 
face. Near the critical temperature, the dynamics of the 
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Figure 3: (Color online) Three different estimators for the 
critical temperature on finite lattices. 
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Figure 5: (Color online) The order-order interface tension 
a(I/) a,t P = 0.318 for different lattice sizes. The dashed line 



is a fit of the data according to Eq. ( 12 1, with fit parameters 
(Too, C2 and C4; the horizontal line indicates the extracted value 
for a no- 
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Figure 4: Entropy density as a function of /3 for L — 56. 



order-order interface aoo is dominated by massless modes 
and its infrared properties are universal (see e.g. pill8|). 
The asymptotic interface tension can then be extracted 
using the ansatz (see e.g. [T^[2U] ) 



<^oo{L) = (To 



C2 



C4 



L2 L4 



(12) 



where the order of the truncation of the expansion in 
1/L^ is determined by the accuracy of the data. As a 
consequence, a reliable extraction of cj requires accurate 
data. From our simulations, we have extracted Fi using 



Eq. (fTl) and then the interface tension fitting the data 
according to Eq. (12 1 truncated to 0(L~*). To reduce 



finite size effects, we included in the fit only points for 
which L^Ja > 6. Since the ansatz ( 12 ) is expected to hold 
only at large distances, our analysis was limited to values 

— 1/2 

of /3 for which Uoo > 3. An example of the quality of 
our data is given in Fig. [5j The values of aoo extracted 
from our fits are plotted in Fig. |6J The relative error on 
this quantity is at most 3 x 10~ , and is invisible on the 
scale of the figure. Near j3c, the behavior of (Too{P) can 
be parametrized as 



Afi) 



,{(i,) + a{fi-p. 



(13) 



Fitting this functional form to our results, we found that 
this provides an excellent description of the data (a fit 
with 9 degrees of freedom has x^/9 = 0.11). We find 
aooiPc) = 0.0249(6) and p = 0.76(4). The quality of the 
fit is shown by the dashed line in Fig. |6] 

At /3c, the order-disorder interface separates a region 
in an ordered state from one in the disordered state. The 
interface tension between an ordered and the disordered 
phase, aod, can be determined by looking at the probabil- 
ity distribution of the energy at the critical temperature. 
In particular, if Pmax is the peak of the histogram when 
the two maxima have equal eight and Pmin is the minimal 
height of the valley between the two peaks [HI [55] , 



1 



2aod{L) = y^log 



P„. 



(14) 



Our data for P{E/N) as a function of E/N are shown 
in Fig. [7| We use a fitting function of the form 



2aod{L) 



(15) 
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Figure 6: (Color online) The infinite volume order-order in- 
terface tension as a function of /3. The fit (dashed line) is 
in excellent agreement with the hypothesis of perfect wetting 
(open triangle). 




Figure 7: (Color online) The probability distribution at the 
critical temperature. The maxima are normalized to one. 



Universality arguments [531 [H] suggest that C3 = 0. Nev- 
ertheless, a correction at this order is due to the fact 
that we use the finite- volume estimator ( 14 ) , as it can 
be shown with a simple saddle point argument [26 . A 



fit of the data according to Eq. (151 performed excluding 
the two smallest volumes yields 2tTod = 0.0252(4), with 
C3 and C4 both compatible with zero. Note that the ob- 
tained value of CTod is fully compatible with the perfect 
wetting relationship 



2CTn 



(16) 



which has been argued in [3S] for the two-dimensional 
Potts model. Our result is a clear indication that perfect 
wetting also holds for the Potts model in three dimen- 
sions. This is shown in Fig. |6] 



V. CONCLUSIONS 

We have proposed a method for determining numer- 
ically free energies of interfaces or topological objects 
when the partition function with given boundary con- 
ditions is required. The range of applicability of our 
method includes not only statistical systems {XY model, 
Heisenberg ferromagnet etc.), but also gauge theories 
(e.g. the 't Hooft loop tension in SU(iV) Yang-Mills). We 
have successfully tested this method on the 3D four-state 
Potts model, for which we have provided a very accu- 
rate determination of the order-order interface below the 
critical temperature. This has enabled us to give clear 
numerical evidence for perfect wetting in this model. 
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